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ABSTRACT 

Aims. In this study, we aim at investigating the properties of the active longitude system of the young solar analogue HD 116956 in 
detail, especially concentrating on determining the rotation period of the spot-generating mechanism with respect to the photometric 
rotation period of the star itself. Because the nonparametric approach, like the Kuiper method, can only give the period of active 
longitudes, we formulate a new method that can determine the parameters the active longitude distribution uniquely. 
Methods. For this purpose, we have developed an analysis method, based on Bayesian statistics using Markov chain Monte Carlo, 
presented in this manuscript. One of the advantages of this method is that an estimate of the active longitude system rotation period, 
as well as the parameters of the shape and location of the active longitudes together with their respective error estimates. This allows 
us to compare the active longitude and mean photospheric rotation periods of the star. 

Results. Our analysis confirms previous results of the object having two stable active longitudes with a phase difference of A(p ss 0.5, 
the other longitude having dominated over the other one during almost the entire span of the time series. Our method gives the rotation 
period of the active longitude system = 7.8412 ± 0.00002d, which is significantly different from the mean photospheric rotation 
period of the star P rot = 7.8 17 ± 0.003d. 

Conclusions. Our analysis indicates that the spot-generating mechanism, manifesting itself as a system of two active longitudes, is 
lagging behind the overall rotation of the star. This behaviour may be interpreted as a nonaxi symmetric dynamo wave propagating in 
the rotational reference frame of the stellar surface. 

Key words. Methods: data analysis, Stars: activity, starspots, individual: HD 116956 



1. Introduction 



In our recent paper dLehtinen et al.ll20TTl) . we presented an anal- 
ysis of 12 years (from late 1998 till mid 2010) of Johnson V 
band photometry from the star HD 1 16956. This star is a chro- 
mospherically active young solar analogue having a spectral type 
of G9 V, and an age of the order of a few 100 Myr (Gaidos ~et al.l 
2000). The spot activity of the star causes rotational modulation 
of brightness, occasionally reaching as high as 0.07 mag. 

Our analysis revealed that the star has stable active longi- 
tudes that have retained their identities over the whole 12 years 
of observations. For the most of this time, one of these active lon- 
gitudes has been dominant, accompanied by a weaker secondary 
that seems to disappear and reappear irregularly. The only ex- 
ception to this behaviour was observed during the first observing 
season in 1998-1999, where we identified a flip-flop, i.e. switch 
of the strongest spot activity from one active longitude to the 
other. The phase separation of the two active longitudes has re- 
mained close to A<p * 0.5. We determined the rotation period 
of this active longitude structure to be Pal = 7.841 6 + 0.001 Id, 
using the formulation of the Kuiper test given by iJetsu & Pelt! 
(fl996h . This value was somewhat different from the weighted 
mean photometric rotation period P lot = 7. 8288d suggesting that 
the active longitudes rotate with a different period than the star 
itself. 

Active longitudes migrating in the rotational reference frame 
of the star ha v e bee n found by the spectroscopic study by 
Lind borg etall (1201 ll) from the RS CVn binary II Peg. In this 
case the active longitudes rotate faster than the orbital period 
of the star in the synchronised binary system. This is in agree- 



ment with t he pho t ometr ic analysis of Berdvugi na & Tuominenl 
(1998) and iJetsul dl996l) . In the photometric studies, migra- 
tory active longitudes were fo und for several other stars, in- 
cluding EI Eri and HR 7275 (IBerdvugina & Tuominenl fl 998) 
and A And and V711 Tau Oetsul 119961) . For t he stars cr 
Gem (IBerdvugina & Tuominen 1998) and FK Com (IJetsu et alJ 
Il993l) . the two periods seem to coincide. Especially in the case 
of FK Com, indications of the migration period b eing variable in 
time i s piling up based on spectroscopic studies (Korh onen et al.l 
120071) . 

In this paper, we examine the behaviour of the active longi- 
tudes of HD 116956 in greater detail. For this purpose, we for- 
mulate a new analysis method based on Bayesian statistics. In 
Sects. |2]and|3] we discuss the background of the statistical anal- 
ysis of active longitudes and the setup of the probability model 
used in the Bayesian analysis. In Sect.|U we formulate the anal- 
ysis algorithm using Markov chain Monte Carlo. This algorithm 
is then applied to the data from HD 116956 in Sect. [5] Finally, 
in Sect. [6] we discuss our results in the light of dynamo theory. 



2. Background 

For the analysis of active longitudes, we need position informa- 
tion of the active areas on the star as a function of time. Minimal 
necessary data consists of central meridian passing times f, of 
the active areas. In the case of photometric time series analysis, 
these are provided by the light curve minimum epochs. If the 
active areas on the star are concentrated on stable active longi- 
tudes, the tj should occur regularly with the active longitude rota- 
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tion period P a i- Longitudinal position information, in the form of 
phases 0, of the inferred active areas, is obtained by folding the 
t, with the period P a i- It is typically well justified to assume the 
number of active longitudes to be between zero and two. This is 
because from the photometry alone it is very difficult to separate 
two simultaneously active areas from each other if their p hase 
separation is smaller than <p as 0.33 (see lLehtinen et al1l201 ll) . 

To determine the value of P a \, the minimum times ?,■ are typ- 
ically analysed with some nonparametric analysis me thod, such 
as the Kuiper test dKuipeJI 19601: iJetsu & Pehl 1 1996h . Such an 
analysis was done by IJetsu! d!996l) for four active stars to find 
active longitudes in their photometry. A more detailed quantita- 
tive analysis can be built on the basis of a preliminary nonpara- 
metric search for P a i by parametrising the full phase distribution 
of (pj. This way we can both improve the precision with which 
P a i is determined, but also quantify the shape and location of 
the active longitudes in the reference frame rotating with the pe- 
riod P a i. The solution is obtained by using the tools of Bayesian 
statistics. 

In this paper, we analyse the light curve mi nimum epochs 
of HP 1 16956 determined with the CPS method dLehtinen et alJ 
1201 lb . The input data consists of all reliable primary and sec- 
ondary minimum epochs t\j and tjj, and their error estimates cr tl . 
and o~ t2i . These epochs are given in Heliocentric Julian Dates. 



3. Setup of the model 

We use Bayesian inference for modelling the active longitudes. 
In other words, we compute the joint posterior probability den- 
sity p(9\y) of the model paramete rs 9 conditional to th e observed 
data y using the Bayes' formula (Gelman et al. 2004) 



wrapped normal distribution (Batschelet 19811) . 
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Unfortunately the wrapped normal distribution cannot be 
written in a closed form. A standard approximation is the von 
Mises distribution (Batschelet 198j]), 



/?(<AI<Ao,k) 



1 



a KCOS (l//-l//o) 



for phase angle ip e [0, 2n]. Here Iq(x) denotes the modified 
Bessel function of the order 0. For the concentration parame- 
ter k — 0, this gives a uniform distribution in phase angle and 
for large k this approaches the wrapped normal distribution with 

We are interested in analysing our observations using phases 
<p rather than phase angles tp. Therefore we rewrite the von Mises 
distribution as 



P (0100, k) = 



1 



cos (2rt<t>-(f>o)) 



(2) 



for phase e [0, 1]. 

One von Mises distribution is adequate for describing one 
active longitude. However, there is often reason to suspect that 
the star being analysed has two active longitudes rather than one. 
These can be modelled by a mixture of two von Mises distribu- 
tions as 



p(<p\0) - m- 



-,/ccos (2^(0-0o.i )) 



/o(*i) 



+ (1 - m)- 



-,/c cos (2jt(0-0 o .2)) 



(3) 



P(0\y) = 



p(0)p(y\Q) 
p(y) ' 



(i) 



The important constituents in the computation of the poste- 
rior density are the prior density p(9) and the likelihood function 
p(y\9). The likelihood function is the actual probability model « 
for the data. It gives the probability density of the observed data p(<p\9) = |~~[ p(<pi\9) 
conditional to given values of the model parameters. The prior 
density p(8) quantifies any existing information about the model 
parameters. This can be, for example, information gathered dur- 
ing previous experiments. The prior can also be formulated so 
that it quantifies the lack of information about the parameter val- 
ues in an objective way. The factor p(y) = J p(8)p(y\9)d9 de- 
pends only on the observed data and acts as the normalisation 
coefficient of the posterior density p(8\y). 

Knowing the full probability distribution for each of the 
model parameters, i.e. the marginal posterior densities, we can 
calculate the parameter estimates as means or medians, and their 
errors as standard deviations or posterior probability intervals of 
the respective distributions. 



where 9 = [m, 0o,i, <Po,2, ft, ""2] are the free parameters of the 
mixture model. To ensure correct normalisation, the mixture pa- 
rameter m has to be restricted to the interval me [0,1]. 

As the data in reality consists of n data points, the actual 
likelihood function is a product of likelihoods for each of these, 



(4) 



It is also possible to use weights for the individual data points, 
in which case the likelihood becomes 



P^<P\e) = \\p(.<Pi\o) w '. 



(5) 



Here the normalised weights w, can be obtained from any rela- 
tive weights vv, by dividing with their mean 



Wi = 



nwi 



(6) 



3.1. Likelihood function 

We are interested in modelling the distribution of phases of the 
light curve minima. Thus, our likelihood function p(y\9) should 
be some circular probability distribution. Let us assume, that 
there is an active longitude on the stellar surface producing ob- 
servable active areas around some phase <po and that there is both 
physical and observational uncertainty in this phase. It is natural 
to think that the phase uncertainty has a Gaussian distribution. 
Using phase angle \jj - 2n<p (ip e [0, 2n]) this is equivalent to the 



Hence, the mean of w, is 1 . This normalisation of the weights is 
important, because the product in Eq. [5] must contain the same 
amount of information from the n available data points as the 
nonweighted product in Eq. [4] 

Lastly, to be able to model the active longitudes from the 
observed light curve minimum epochs, we need the active longi- 
tude period P a i (hereafter called P when used to denote the free 
model parameter) to convert these time points f, into phases <pj. 
This is done by 



•i(P) = FRAC 



(7) 
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where FRAC[x] removes the integer part of x. The period P need 
not to be given, but can be considered as the sixth free model 
parameter together with the other free parameters 9 of Eq. [3] 
The period P is included in the likelihood function implicitly via 
Eq.El 

3.2. Parameter priors 

In order to compute the posterior density p(9\y), we need priors 
for the model parameters 6. We choose noninformative priors, 
i.e. priors that quantify our lack of knowledge about the values 
of these free parameters. 

For the active longitude phases 0o,i and 0o,2, the most nat- 
ural choice for a noninformative prior is a uniform prior in the 
interval [0, 1]. This states that a priori we consider that the active 
longitudes may lie at any phase with equal probability. Similarly, 
the mixture parameter m is bound to the interval [0, 1] and is nat- 
urally given the same uniform prior. Note, however, that the pri- 
ors p(0o,i) and p{<poj) are only formally bound between and 1. 
In actual modelling, we allow 0o,i and 0y,2 to get values below 
as well as above 1 . The sole purpose of this approach is to take 
care that the wings of the marginal posterior densities of 0o,i and 
0o,2 would not be cut at the phases <p — or 1 . 

For the concentration parameters k\ and k 2 , we 
choose the improper, i.e. n on-normalisable, Jeffreys' prior 
dDobigeon & Tourneretll2007l) 

p(k) = -, (8) 

K 

for k > 0. This is constructed specially to reflect vague infor- 
mation about the value of k. The same prior works also for the 
period dGregorv & Loredoll992l) . i.e. 

P(P) = ^, (9) 

for P > 0. The fact that this prior is noninformative about P 
can be seen by a transformation of variables / = P l from the 
period domain to the frequency domain. Under this change, the 
prior given by Eq. [9]retains its shape, just as a noninformative 
prior should. 

The joint prior of all the model parameters is just the product 
of all the individual parameter priors, 

p(8) = p(P)p(m)p((f> OA )p((p oa )p(Ki)p(K2). 

4. Modelling algorithm 

To estimate the joint posterior density p{6\y) of the model param- 
eters, we need to use a Monte Carlo approach. In other words, 
we simulate random numbers that are distributed according to 
p(0\y). This is possible by usi ng the Markov ch ain Monte Carlo 
(hereafter MCMC) approach (iGilks et al.lll996b . 

4.1. Preliminary estimation of P 

To ensure that the MCMC algorithm indeed converges to the 
correct parameter values, it is necessary to make a preliminary 
search first for the posterior mode within broad search intervals. 
Because our model (Eqs.[3}|7]l is rather complex, this search is 
done in several stages. 

The first model parameter to be searched for is the active lon- 
gitude period P. For this search w e use the Kuiper periodogram 
(lKuiperlll960l:lJetsu & Pelll[T996l) . constructed by computing the 



value of the Kuiper test statistic V„(f) for a discrete grid of fre- 
quency values / = P . 

Computing the Kuiper test statistic consists of comparing the 
empirical cumulative distribution function F„(4>) of the phases 0, 
to some theoretical cumulative distribution function F(<p). As we 
are interested in finding any deviations from uniform distribution 
of <f>i, the theoretical cumulative distribution function in our case 
is F(<p) = <p. The test statistic is 

V„(f) = D + + D r 

where D + and D denote the maximum values of F„{<p) - F((f>) 
and F(<p) - F„((p), respectively. 

As stated above, the Kuiper periodogram is computed for a 
fixed set of frequency values. This is most naturally done with a 
grid centered around the initial guess for the frequency /o = P^ 1 . 
The tested frequencies are within the range 

(1 -q)fi> </<(l+?)/o, 

and are separated by evenly spaced steps of 

Af=[(t n -h)OFACT i 

from each other, where t\ and t„ are the first and the last time 
point, respectively, and OF AC is the over filling factor. In this 
paper, we have used q = 0.05 and OF AC = 10. The period 
value P - f~ l which maximises the periodogram, i.e. gives the 
largest value of V„(f), is used as the preliminary estimate for P. 

4.2. Preliminary estimation ofm, o ,i. <t>o,2, ki and k 2 

With a preliminary estimate for P, we can start searching for 
estimates for the rest of the model parameters. This can be con- 
sistently done with the full MCMC strategy, by finding the ap- 
proximate mode of the posterior density (Eq. Q3 conditional to 
the value obtained for P. 

To speed up the process of finding the approximate posterior 
mode, it is useful to divide the parameters into two groups. We 
first search for the modes of m, 0o,i and 0o,2, each within the in- 
terval [0,1]. This is done by computing the value of the posterior 
density with proposed values of m, 0o,i and 0o,2, using prelim- 
inary estimate for P and some reasonable initial guesses for k\ 
and k 2 , and taking steps towards higher posterior density. Once 
reasonable estimates for m, tpo,\ and 0o,2 have been obtained, we 
search for the modes of k\ and k 2 . This is done by computing 
the posterior density for a grid of k\ and k 2 conditional to the 
preliminary estimates of the rest of the parameters. In this paper, 
we do the initial search between K m i n = 0.5 and K max = 50.0 and 
with a grid spacing of Ak = 0. 1 . 

4.3. Posterior simulation 

After crude estimates of the model parameters have been ob- 
tained, the full po sterior density p (9\y) is modelled with a 
MCMC algorithm (Gil ks et alj|1996h . This algorithm generates 
random values for the model parameters as Markov chains. In 
other words, we simulate «mcmc values for all of the model pa- 
rameters 8, so that the values on each simulation round, 6,, de- 
pend on their values on the previous round, . When the algo- 
rithm is executed correctly, these chains converge to the marginal 
posterior densities of the individual model parameters. The pre- 
liminary parameter estimates are used as the starting values for 
the Markov chains. 

Because of the complex nonlinearity of the active longitude 
model (Eqs. [3]-[7j, we split the MCMC algorithm into several 
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stages. During these, a limited number of free parameter values 
are updated. 

The first parameter to be updated is the mixture parameter 
m. For this, we compute the value of the posterior density at 
each minimum phase with the parameter values of the last sim- 
ulation round and assuming likelihoods containing only one ac- 
tive longitude, i.e. either = p(0,j0o,i>*i) or p{<pi\2) = 
/>(0;l0o,2, K2). Using these, we compute the probability of 0,- be- 
longing to the active longitude 1 , 



Zi = 



m t _i/?(0,-|l) + (1 - ot ; _i)p(0,|2)' 



(10) 



where m,_i denotes the value of the mixture parameter on the 
previous simulation round. The new value for m is now given by 
the weighted mean of Zi for all </>,-, 



L i= i win 



(ii) 



This procedur e for simulating m is based on the strategy pro- 
posed bv lGelman et al.l (|2004) for computing mixture models. 

The computation of the probabilities z, also allows us to se- 
lect the active longitude that the minimum phases <p, belong to. 
The division is done at each simulation round and is probabilis- 
tic, i.e. each 0, is assigned to the active longitude 1 with the 
probability z,. Those minimum phases that are not assigned to 
the first active longitude are automatically assigned to the second 
active longitude. This division of the minimum phases makes it 
possible to model both of the active longitudes separately. 

After a new value for m has been simulated, we update 
the values of the period P and the phases of both of the ac- 
tive longitudes, 0o,i, 00,2, using the bimodal likelihood function 
(Eq. [3]) with the new value for m and old values for k\ and K2- 
The simul ation is done wit h a standard Metropolis-Hastings al- 
gorithm dGilks et al.| [T996), where new proposed values 8, for 
m, 0o,i and 0o,2 are drawn from a Gaussian proposition den- 
sity q(0,) = N{6,\6,-\,<T 2 qg ) centered around the parameter value 
6,-i of the previous simulation round and having a different ad- 
justable variance cr 2 g for the different parameters. 

Last, we update the values of the concentration parameters 
K\ and K2- Since obtaining good mixing for the Markov chains of 
these parameters can be tricky, we do the updating separately for 
both of these parameters using a unimodal likelihood function 
(Eq.|2]i and only the minimum phases labelled as those belonging 
to the active longitude of the currently updated concentration 
parameter. 

We use a gamma distribution as the proposition density for 
the conce ntration parameters, q(K,) = Gamma(K,\a,/3), as sug- 
gested by Dobigeon & Tourneretl(l2007l) . The parameters for this 
distribution are set to be a = k 2 /<t 2 K and j3 = k/cr 2 K , where k de- 
notes the preliminary estimates for the concentration parameters 
and cr 2 K is an adjustable variance for the proposition density. 

It usually takes some amount of simulation rounds before 
the Markov chains converge to the marginal posterior densities 
of the model parameters. To eliminate this burn-in phase in the 
beginning of the simulation, it is necessary to discard the first 
"bum-in simulation rounds from the chains. In short simulations, 
it may be necessary to set nbum-in to be 50 % of the total length 
of the simulation. In longer simulations, 10 % of the total length 
is often sufficient. 

The marginal posterior densities of the model parameters 
might not be very illustrative as themselves. Therefore we sum- 
marise the results of the MCMC simulation by computing esti- 
mates for the parameter values and their errors. If the parameters 



Table 1. Parameter and error estimates for global analysis of the 
active longitudes. The second column gives the means and stan- 
dard deviations of the Markov chain of P, <pi and <p2- The third 
and fourth columns give the medians and 5th and 95th centilles 
of the Markov chains of all of the parameters. The active lon- 
gitude rotation period P is given in days while the rest of the 
parameters are dimensionless. 





mean ± std 


median 


[5%, 95%] 


Pal [d] 


7.8412 ± 0.0002 


7.8412 


[7.8409, 7.8414] 


111 




0.78 


[0.75, 0.80] 


0. 


0.0 ± 0.007 


0.0 


[-0.011, 0.011] 


4>2 


0.479 ± 0.019 


0.479 


[0.444, 0.511] 


Kl 




4.88 


[4.28, 5.62] 


«2 




1.21 


[0.75, 1.86] 



turn out to have nearly Gaussian distributions, it is natural to esti- 
mate their values with the means and their errors with their stan- 
dard deviations of the respective Markov chains. For more com- 
plicated distributions, it is better to give the median and some 
quantiles q and 1 - q of the Markov chains as the parameter and 
error estimates. 



5. Analysis of the active longitudes 

We use the MCMC algorithm described above to analyse the ac- 
tive longitudes of HD 116956. The original observations were 
made between HJD = 2451172 (24th of December, 1998) and 
HJD = 2455342 (25th of May, 2010) with the T3 0.4 m auto- 
matic photoelectric telescope (APT) at the Fairborn Observatory 
in Arizona. They were divided into 12 continuous observing sea- 
sons, each spanning roughly the first half of the ye ar. We labelled 
these observing seasons as segments SEG 1—12 dLehtinen et al.l 
201 1). The CPS method provided estimates for both the primary 
and secondary minimum epochs, t \ j and t2j. However, in our 
analysis of this paper we pool these two classes together as f,, 
where the number of time points is n — 645. This pooling is done 
because, even if the star has two stable active longitudes, both the 
primary and secondary light curve minima can very well appear 
on both of the active longitudes. This is possible if the levels of 
spot activity on the two active longitudes are close to each other 
and either one of them may be stronger at any given moment. 
Another possible alternative is a flip-flop, where the main spot 
activity physically shifts from one active longitude to the other. 

For the analysis in this paper, we use weights. In other words, 
we associate a weight w, for each data point t, and use Eq.[5]to 
combine the likelihoods of these data points. A typical definition 
for relative weights is the inverse square of error, i.e. w, = <r~ 2 . 
This is also what we use in the CPS. It may, however, occur that 
the weights obtained this way have a large dispersion and some 
of them have values greatly larger than others. In this case, the 
burden of carrying information from the data to the posterior 
density is strongly concentrated on just a few data points. This 
can in turn lead to an instability in the posterior simulation. To 
even out the pathologically defined weights, we must use a more 
robust definition for them, e.g. w, - <rj l . We use this definition 
for the weights throughout the paper. 

5.1. Global and seasonal analyses 

We are interested in both the long term overall behaviour, as 
well as any seasonal changes of the active longitudes. Thus, 
we present both a global analysis, including all the available 
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Fig. 1. The simulated parameter values of the MCMC run for the 
global analysis. Presented are scatter plots between each pair of 
model parameters showing the correlations between the param- 
eters. 



minimum epochs, and seasonal analyses that include minimum 
epochs only from one observing segment at a time. For each 
case, we ran the MCMC algorithm with hmcmc = 

10000 rounds 

and discarded the first Hbum-in = 1000 as the burn-in phase. 

An idea of how the MCMC algorithm works in practice, is 
visualised in Fig. Q] It shows the Markov chains of the model 
parameters 9 in the global analysis as scatter plots between pairs 
of free parameters. As the burn-in phase has been discarded, 
these scatter plots should reflect the shape of the multidimen- 
sional joint posterior density. In an ideal case, each panel in Fig. 
Q]should approximate a bivariate Gaussian with vanishing corre- 
lation between the parameters. In reality, there are some param- 
eter pairs, such as m and K2, which show relatively strong corre- 
lation or anticorrelation. This shows that the parametrisation of 
the problem is not the best possible. It is, however, not clear what 
kind of a reparametrisation could remove this correlation. In any 
case, the present more or less "physical" parametrisation seems 
also to yield good mixing of the Markov chains, i.e. the joint 
posterior density is uniformly sampled. In addition to significant 
correlation, some panels of Fig. [T]show slightly irregular shapes. 
Thus, the joint posterior density has a more complex form than 
what could have been expected a priori. 

Parameter estimates of the global analysis are given in Table 
Q] Because the marginal posterior densities of m, k\ and K2 are 
asymmetric, i.e. not Gaussian, we estimate the parameter val- 
ues with the medians of the Markov chains of each parameter. 
For error estimates, we give the [5%, 95%] posterior intervals. 
The marginal posterior densities of <p\, <fa and P, now denoted 
with P a i = P being regarded as the rotation period of the ac- 
tive longitude system, have more closely Gaussian forms, and 
for them, we also give the mean and standard deviation of the 
Markov chains. The two estimates based on the mean or median 
are in good agreement with each other. 




Fig. 2. The distribution of all light curve minima folded with the 
ephemeris HJD al = 2451177.6135 + 7.48 12E (Eq.QJi and pre- 
sented as a histogram between the phase interval -0.5 < <p < 0.5. 
On top of the observational distribution, we show the active lon- 
gitude model distribution using the estimated model parameters 
from Table Q] 

Based on the global analysis, we calculated a new linear 
ephemeris for the central meridian passing epochs of the primary 
active longitude (labelled with 1), 

HJDai = 2451177.6135 + 7.4812£, (12) 

where E is the number of rotations. The model distribution (Eqs. 
[3] and [7]i using the estimated parameter values is compared with 
the observational distribution of the light curve minima in Fig. 
[2] folded with the period P d \ = 7.4812 d. As suggested by the 
value m = 0.78, the primary active longitude is clearly the dom- 
inant structure on the star with the secondary active longitude 
only barely visible. The phase difference between the two active 
longitudes is quite close to A(p = 0.5. The model fit to the ob- 
servational phase distribution is reasonably good, although not 
perfect. The differences may be due to seasonal changes in the 
active longitude phases. 

The seasonal changes of the active longitudes were investi- 
gated by setting the period P - P a \ - 7.8412 d and not letting 
this value vary in the MCMC runs of the individual segments. 
Parameter estimates for each of the 12 segments are presented 
in Table [2] The phases are consistent with the ephemeris of Eq. 
[T2l and therefore the primary active longitude phase <p\ has val- 
ues near both and 1. Table [2] also includes the phase difference 
A(p of the primary active longitude and the following secondary 
active longitude, when the two have been present. 

The seasonal active longitude phases, <p\ and <p2, are plot- 
ted in Fig. [3] along with the phases from the global analysis. It 
is immediately clear, that the active longitudes have undergone 
significant shifts. This is in good agre ement with our previous 
results (Fig. 10 in lLehtinen et al.ll201 lb . In the present analysis, 
also the phase difference A<p between the active longitudes seems 
to vary significantly. This may, however, be caused by the small 
number of secondary light curve minima available in many of 
the segments, giving a less stable simulation of the secondary 
active longitude phase. 

In the segments 5, 6, 8, 10 and 12, the MCMC algorithm 
yields mixture parameter values of m = 1, meaning that only 
the primary active longitude is present. The algorithm does not 
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Table 2. Parameter and error estimates for local segment by seg- 
ment analyses of the active longitudes presented in the same 
form as in Table [TJ The phase difference A<p between the pri- 
mary and following secondary active longitude is also given, 
when both of the longitudes are present. The active longitude 
period of the global analysis, P d \ = 7.8412 d, has been used for 
all of the segments. The year during which the largest part of the 
segment data was observed is also given. 
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Fig. 3. Active longitude phases folded with the ephemeris 
HJD al = 2451 177.6135 + 7.4812£ for primary (black symbols) 
and secondary (grey symbols) active longitudes. Individual sym- 
bols denote seasonal and horisontal lines global estimates. The 
lcr errors are denoted with error bars for the seasonal and with 
dashed lines for the global estimates. 

Table 3. Seasonal active longitude periods with lcr errors. 



SEG 


Pal M 




1 


7.8871 


± 0.0220 


2 


7.7031 


±0.0136 


3 


7.8067 


± 0.0077 


4 


7.8952 


±0.0115 


5 


7.8147 


±0.0134 


6 


7.8651 


± 0.0080 


7 


7.8752 


± 0.0097 


8 


7.9128 


± 0.0046 


9 


7.8555 


± 0.0036 


10 


7.8193 


± 0.0047 


11 


7.7695 


±0.0131 


12 


7.8300 


± 0.0223 



determine the number of active longitudes present in the data, 
but uses the bimodal model of Eq.[3]for all data. In spite of this, it 
was able to detect the absence of the secondary active longitude 
in these five segments. 

Lastly, it is worth noting that the concentration parameters ki 
and K2 have very different values in individual segments and typ- 
ically also wide [5%, 95%] posterior intervals within each seg- 
ment. This demonstrates the fact that these parameters are the 
least stable ones of our model. The typical problem of border- 
line data is often just to get decent mixing for the Markov chains 
of Ki and K2- This is, however, not a major concern, since the ex- 
act values of k\ and K2 are of little interest. Our main goal is to 
determine the values of P, <j>\ and 02- 

5.2. Comparison with the mean rotation period 

The global constant value of P a \ was used for modelling each of 
the segment in the previous section. We perform another analy- 
sis, where the period of each segment is also an additional sixth 
free parameter, just as in the case of the previous global analysis. 
The active longitude phases obtained this way are not consistent 
with each other, but this approach can give information of the 
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Table 4. Global and seasonal weighted mean photometric peri- 
ods with lcr errors. 



SEG 


P mt [d] 




Global 


7.817 ± 


0.003 
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7.890 ± 


0.014 
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7.698 ± 


0.014 
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7.808 ± 


0.010 
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7.803 ± 


0.021 
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7.801 + 


0.009 
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7.855 ± 


0.010 
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7.830 ± 


0.011 
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7.836 ± 


0.008 
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7.817 ± 


0.005 
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7.807 ± 


0.010 
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7.839 ± 


0.017 


12 


7.659 ± 


0.044 



changes in the active longitude rotation period. The resulting pe- 
riod estimates presented as means and standard deviations of the 
Markov chains of P are given in Table [3] 

We want to compare the global and seasonal periods of the 
active longitudes computed in Sect. 15. ll with the mean period of 
the light curve of HD 1 16956. The idea is that the photometric 
periods obtained with the CPS correspond to the rotation of the 
active areas on the stellar surface. On the other hand, active lon- 
gitudes that stay stable for a number of years would rather seem 
to represent the rotation of the underlying global geometry of the 
magnetic field producing the active areas. 

The two kinds of rotation periods need not to coincide. If 
the rotation period of the magnetic structure significantly differs 
from the surface rotation period, we would expect to observe the 
active areas forming at the active longitudes, but then start to 
drift a way. 

In iLehtinen et alj d201 lb we calculated the weighted mean 
light curve period of HD 116956 from the results of the CPS 
analysis. To be able to compare the mean period with the ac- 
tive longitude rotation period, we also need an error estimate for 
the former. Such a weighted mean period and its error estimate 
can be computed in the Bayesian paradigm consistently with the 
active longitude analysis presented above. We thus need to for- 
mulate the posterior density of the mean photometric period P mt . 

Each period estimate P, from the CPS can be assumed to 
have a Gaussian error distribution with the variance cp- . This 
corresponds to the likelihood 

piP^o-l) = N(Pi\P mt ,o\) = : 1 er^-^l^l,. 

crp. y2n 

Combined for all the period estimates, this becomes 

n 

p(P\P m[ ,cr 2 p ) = Y\X(Pi\Pmucrj) 

i=l 

(27r)"/2 1 1 L p . J 
By using the prior of Eq.[9]for P lot , the posterior density becomes 

p(P mt \P,oi) cc p(P I0t )p(P\P wU oi) cc P ro 1 t e Z;Ll [" (P '°'^ )2/2 ^].(13) 

As each P, here has a different crj; , the estimate of P rol using Eq. 
[T3lis indeed weighted. 

Estimates for both the global and seasonal weighted photo- 
metric mean periods were obtained by sampling the posterior 
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Fig. 4. Variations of the seasonal estimates of the active longi- 
tude rotation period P a i (black squares) and the seasonal photo- 
spheric rotation period P lot (grey triangles). The long-term active 
longitude period from the global analysis and the weighted mean 
photospheric rotation period are presented with black and grey 
horisontal lines, respectively. The error bars of the seasonal es- 
timates and dashed lines of the global estimates denote lcr error 
limits. 

density of Eq. [13] with the Metropolis-Hastings algorithm. For 
the sampling, we used a Gaussian proposition density q{P mt ,r) = 
N{P m u\Pmi,t-\,cr 2 qV ) with cr fyPioi = 0.01 d. For each estimate, 
we ran the MCMC algorithm with «mcmc = 4000 rounds, first 
«bum-in = 2000 of which were discarded as the burn-in phase. 
The period estimates presented in Table|4]are means of the sim- 
ulated Markov chains and the respective standard deviations are 
assigned as their error estimates. 

The results are compared graphically with the correspond- 
ing active longitude rotation periods in Fig. |4] Both the photo- 
metric and active longitude rotation periods seem to have varied 
quite much during the years. For the first years of observations, 
the two types of seasonal periods seem to have coincided rather 
closely, but this behaviour is lost in the more recent observations 
and the two periods show no mutual correlation. This supports 
the view that these two periods trace different physical mecha- 
nisms on the star. 

When comparing the global rotation periods of the active 
longitudes, P d \ = 7.8412 + 0.0002 d, and the photospheric ac- 
tive areas, P mt = 7.817 + 0.003 d, the difference between the 
two is clear. The long-term mean photospheric rotation period is 
significantly shorter than the active longitude rotation period. In 
the seasonal period estimates, this difference is not that clear be- 
cause of their large variability. Nevertheless, the effect is visible 
as all except one seasonal mean photometric periods are shorter 
than their corresponding active longitude periods or equal within 
the error limits. Thus it seems, that the active longitudes of HD 
1 16956 rotate more slowly than its photosphere. 

6. Discussion in the light of dynamo theory 

It is generally believed that the magnetic fields in the active rapid 
rotators are due to a dynamo mechanism working almost solely 
on the inductive action arising from convective turbulence. The 
other basic ingredient arising from differential rotation, i.e. the 
Q-effect, is expected to be negligibly weak in the rapid rota- 
tion regime due to the r elative differential rotation A Q/Q be- 
ing strongly reduced (e.g. Kitchatinov & Rudiger 1999). Such a 
dynamo is called as an a 2 -dynamo, the a-effect describing the 
collective inductive action of turbulence. 

In the regime of slower rotation, i.e. in the case of the Sun, 
the O-effect driven by the angular velocity gradients is strong 



8 



J. Lehtinen et al,: Analysing the active longitudes of the young solar analogue HD 1 16956 using Bayesian statistics 



compared to the a-effect. For this type of an o-Q dynamo, oscil- 
lating axisymmetric, i.e. field not changing over the azimuthal 
(longitudinal) coordinate, configurations are typical. The solar 
nearly dipole field showing a magnetic cycle of roughly 22- 
years, with a related latitudinal dynamo wave forming the butter- 
fly diagram, is thought to be a represent ative case of t he action 
of an aQ-dynamo, see e.g. the revi ew bvlOssendriiverl (2003). 

According to both li near (e.g. | Kr ause & Radler 198CJ) and 
nonlinear solutions (e.g. iMoss et al.l 1 19951) of the q 2 -dvnamo 
equations, the nonaxisymmetric modes become more easily ex- 
cited in the rapid rotation regime. The m = 1 mode, represent- 
ing an azimuthally varying field changing sign once over the 
full longitude span, is commonly the preferred field configu- 
ration. This is not surprising, as the largest scale mode is al- 
ways the most resistive over diffusive effects. The nonaxisym- 
metric modes turn out to be waves migrating in azimuthal direc- 
tion, not necessarily ha ving the rotation period of the star (e.g. 
iKrause & Radler| [l980). Both slower and faster dynamo waves 
can occur, depending for example on the profile and properties of 
the turbulent transport coefficients. The faster waves with dipole 
symmetry (SI) were observed to be preferr ed in linear models 
with simple profiles (Krause & Radler 1980), and slower waves 
with quadrupolar symmetry (Al) in mor e complicated nonlinea r 
models solving also for the dynamics (Tuomin en et al.l 12002). 
From the viewpoint of dynamo theory, therefore, two migratory 
active longitudes are an expected result. 

Another peculiar property of the dynamo solutions in this 
regime include their steadiness, i.e. no oscill atory solutions ca n 
be found in the very simplest settings (e.g. IMoss et al.l [T99 ll) . 
Therefore, time-dependent phenomena, such as the flip-flop, are 
difficult to explain in the framework of a simple a 2 dynamo 
mechanism. The usual remedies are to include some differential 
rotation, t hat makes time-dependent solutions more easy to ex- 
cite (e.g. Moss 2005), or to use more compli cated profiles f or the 
turbulent transport coefficients (e.g. Riidiger et al. 2003), nor- 
mally working in the same direction. As the observations cannot 
yet tell whether a polarity change is related to the flip-flop cy- 
cle, it is not yet possible to pin down more exactly which type 
of a dynamo (the possibilities being a pure a 2 or a differential 
rotation aided a 2 Q. dynamo), not to mention which mode of its 
solutions, is responsible for the observed magnetic fields in rapid 
rotators. 



fits t o the phases of the active regions in the reference frame 
(e.g.[Berdvugi na~& Tuominenll998 : lKorhonen et alJ20 07). This 
often involves some subjectivity in choosing which light curve 
minima correspond to which active longitude, as well as making 
use of phases well outside of the range € [0, 1]. 

Our analysis of H D 116956 confirms our previous results 
dLehtinen et al J 120111) . as we find two stable active longitudes 
with a mutual phase difference of A0 * 0.5. One of these 
active longitudes has remained nearly constantly stronger and 
dominated the overall distribution of active areas on the stel- 
lar surface. Our new estimate for the ephemeris of the central 
meridian passing of the primary active longitude is HJD a i = 
2451177.6135 + 7.4812£ (Eq.IHli. 

We compared the long-term active longitude period P d \ = 
7.8412 ± 0.0002 d with the long-term mean photometric rotation 
period P mt = 7.8 17±0.003 d. The estimates for these two periods 
are markedly different in such a way that the active longitudes 
lag behind the photospheric rotation. This behaviour may be in- 
terpreted as a nonaxisymmetric dynamo mode manifesting itself 
as a longitudinal dynamo wave propagating in the rotational ref- 
erence frame of the stellar surface. 

It should be noted that the comparison between the periods 
P a i and P m assumes that the mean photometric rotation period 
P rot actually measures the stellar rotation period. In the case of 
weak differential rotation this clearly is so; also theoretical mod- 
els indicate weak relative differential rotation for rapid rotators. 
In addition, for synchronously rotating binary systems, which is 
the case for many RS CVn stars, we can associate the orbital 
period P or b with the stellar rotation period. However, even in 
the case of significant differential rotation, as our previous study 
suggests for HD 1 16956, we are confident in using P mt as the ef- 
fective measure for the stellar rotation. This is because it tracks 
the mean rotation period in the latitude zones where the spot ac- 
tivity actually takes place, and where it should be compared with 
the active longitude rotation period P a \. 
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7. Conclusions 

We have analysed the active longitudes of the young solar ana- 
logue HD 116956 in depth. For this purpose, we have formu- 
lated a new analysis algorithm based on Bayesian statistics. This 
method uses a mixture of two von Mises distributions to model 
the phase distribution of the active areas manifesting themselves 
as photometric light curve minima. 

Our analysis algorithm allows us to characterise the active 
longitudes in a systematic way. Both the rotation period P a i 
and the shape and location parameters of the active longitude 
distribution are parametrised in the model and no subjectivity 
is involved in their determination. This represents an advance 
from previously used methods. Typically the rotation period of 
the active long itudes has been determined with nonparametric 
methods (e.g. Jetsu 1996). The nonparametric methods can only 
determine the period, i.e. they do not determine any parame- 
ters that characterize the shape or location of the active lon- 
gitude distribution. Alternatively, the migration rate of the ac- 
tive longitudes with respect to the stellar rotational reference 
frame can be studied. This has been done by computing linear 
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